GAMMs at maximum constriction
splines.it$c2phonation.ord <- as.ordered(splines.it$c2phonation)
splines.it$c2place.ord <- as.ordered(splines.it$c2place)
splines.it$vowel.ord <- as.ordered(splines.it$vowel)
contrasts(splines.it$c2phonation.ord) <- "contr.treatment"
contrasts(splines.it$c2place.ord) <- "contr.treatment"
contrasts(splines.it$vowel.ord) <- "contr.treatment"
max.it <- filter(splines.it, label %in% c("max_TT", "max_TD")) %>%
na.omit()
max.it.gam <- bam(
Y ~
c2phonation.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr"),
data = max.it,
method = "ML"
)
summary(max.it.gam)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord,
## bs = "cr")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -4.8549 0.1839 -26.406 < 2e-16 ***
## c2phonation.ordvoiceless 1.3251 0.2592 5.113 3.27e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 8.271 8.772 774.019 < 2e-16 ***
## s(X):c2phonation.ordvoiceless 1.944 2.464 4.707 0.00501 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.676 Deviance explained = 67.6%
## -ML = 21185 Scale est. = 95.775 n = 5721
max.it.gam.null <- bam(
Y ~
s(X, bs = "cr"),
data = max.it,
method = "ML"
)
compareML(max.it.gam, max.it.gam.null)
## max.it.gam: Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord,
## bs = "cr")
##
## max.it.gam.null: Y ~ s(X, bs = "cr")
##
## Chi-square test of ML scores
## -----
## Model Score Edf Chisq Df p.value Sig.
## 1 max.it.gam.null 21203.08 3
## 2 max.it.gam 21185.13 6 17.943 3.000 7.912e-08 ***
##
## AIC difference: -31.72, model max.it.gam has lower AIC.
plot_smooth(max.it.gam, view = "X",
plot_all= "c2phonation.ord", rug = FALSE)
## Summary:
## * c2phonation.ord : factor; set to the value(s): voiced, voiceless.
## * X : numeric predictor; with 30 values ranging from -44.040500 to 61.699900.

plot_diff(max.it.gam, view = "X",
comp = list(c2phonation.ord = c("voiceless", "voiced")))
## Summary:
## * X : numeric predictor; with 100 values ranging from -44.040500 to 61.699900.

##
## X window(s) of significant difference(s):
## -44.040500 - -2.385191
acf_plot(resid(max.it.gam), split_by=list(max.it$rec.date))

max.it.gam.2 <- bam(
Y ~
c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr") +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord), bs = "cr",
data = max.it,
method = "fREML"
)
summary(max.it.gam.2)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord,
## bs = "cr") + s(X, by = vowel.ord)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.3151 0.2486 -25.399 < 2e-16 ***
## c2phonation.ordvoiceless 1.3276 0.2281 5.820 6.21e-09 ***
## c2place.ordvelar 2.3702 0.2340 10.128 < 2e-16 ***
## vowel.ordo 1.8548 0.2777 6.679 2.64e-11 ***
## vowel.ordu -0.2606 0.3050 -0.854 0.393
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 7.728 8.371 414.093 < 2e-16 ***
## s(X):c2phonation.ordvoiceless 2.765 3.469 8.285 7.17e-06 ***
## s(X):c2place.ordvelar 8.112 8.642 102.386 < 2e-16 ***
## s(X):vowel.ordo 5.395 6.592 15.107 < 2e-16 ***
## s(X):vowel.ordu 8.666 8.956 30.062 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.75 Deviance explained = 75.2%
## fREML = 20473 Scale est. = 73.809 n = 5721
max.it.gam.null <- bam(
Y ~
c2phonation.ord +
c2place.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr") +
s(X, by = c2place.ord, bs = "cr"),
data = max.it,
method = "fREML"
)
compareML(max.it.gam.2, max.it.gam.null)
## max.it.gam.2: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord,
## bs = "cr") + s(X, by = vowel.ord)
##
## max.it.gam.null: Y ~ c2phonation.ord + c2place.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord,
## bs = "cr") + s(X, by = c2place.ord, bs = "cr")
##
## Chi-square test of fREML scores
## -----
## Model Score Edf Chisq Df p.value Sig.
## 1 max.it.gam.null 20751.11 9
## 2 max.it.gam.2 20473.03 15 278.072 6.000 < 2e-16 ***
##
## AIC difference: -588.81, model max.it.gam.2 has lower AIC.
plot_smooth(max.it.gam.2, view = "X",
plot_all= "c2phonation.ord", rug = FALSE)
## Summary:
## * c2phonation.ord : factor; set to the value(s): voiced, voiceless.
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 30 values ranging from -44.040500 to 61.699900.

plot_diff(max.it.gam.2, view = "X",
comp = list(c2phonation.ord = c("voiceless", "voiced")))
## Summary:
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 100 values ranging from -44.040500 to 61.699900.

##
## X window(s) of significant difference(s):
## -44.040500 - -8.793700
max.it.gamm <- bam(
Y ~
c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr") +
s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord), bs = "cr",
data = max.it,
method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs", xt = "cr",
## m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + s(X,
## by = vowel.ord)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.3692 1.6283 -3.911 9.29e-05 ***
## c2phonation.ordvoiceless 1.4895 0.4027 3.699 0.000218 ***
## c2place.ordvelar 2.3929 0.4242 5.641 1.77e-08 ***
## vowel.ordo 1.7681 0.5068 3.489 0.000489 ***
## vowel.ordu -0.2405 0.5289 -0.455 0.649397
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 7.338 8.048 20.434 < 2e-16 ***
## s(X):c2phonation.ordvoiceless 2.520 3.103 2.951 0.0297 *
## s(X,rec.date) 261.634 895.000 0.783 < 2e-16 ***
## s(X,speaker) 4.209 8.000 12.604 < 2e-16 ***
## s(X):c2place.ordvelar 8.243 8.709 86.897 < 2e-16 ***
## s(X):vowel.ordo 5.596 6.768 7.128 2.89e-08 ***
## s(X):vowel.ordu 8.675 8.949 23.533 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.784 Deviance explained = 79.6%
## fREML = 20286 Scale est. = 63.649 n = 5721
max.it.gam.null <- bam(
Y ~
c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr") +
s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord), bs = "cr",
data = max.it,
method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
compareML(max.it.gamm, max.it.gam.null)
## max.it.gamm: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, speaker, bs = "fs", xt = "cr",
## m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") + s(X,
## by = vowel.ord)
##
## max.it.gam.null: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") +
## s(X, by = vowel.ord)
##
## Chi-square test of fREML scores
## -----
## Model Score Edf Chisq Df p.value Sig.
## 1 max.it.gam.null 20320.94 18
## 2 max.it.gamm 20286.36 21 34.583 3.000 6.441e-15 ***
##
## AIC difference: -69.31, model max.it.gamm has lower AIC.
plot_smooth(max.it.gamm, view = "X",
plot_all= "c2phonation.ord", rug = FALSE)
## Summary:
## * c2phonation.ord : factor; set to the value(s): voiced, voiceless.
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 30 values ranging from -44.040500 to 61.699900.
## * rec.date : factor; set to the value(s): 12/12/2016 14:51:50.
## * speaker : factor; set to the value(s): sdc02.

plot_diff(max.it.gamm, view = "X",
comp = list(c2phonation.ord = c("voiceless", "voiced")))
## Summary:
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 100 values ranging from -44.040500 to 61.699900.
## * rec.date : factor; set to the value(s): 12/12/2016 14:51:50.
## * speaker : factor; set to the value(s): sdc02.

##
## X window(s) of significant difference(s):
## -44.040500 - -9.861785
acf_plot(resid(max.it.gamm), split_by=list(max.it$rec.date))

Autoregressive correction
max.it <- arrange(max.it, rec.date, fan) %>%
mutate(rec.date = as.character(rec.date))
previous <- ""
for (i in 1:nrow(max.it)) {
if (max.it$rec.date[i] != previous) {
max.it$start.event[i] <- TRUE
previous <- max.it$rec.date[i]
} else {
max.it$start.event[i] <- FALSE
previous <- max.it$rec.date[i]
}
}
r1 <- start_value_rho(max.it.gam.2)
max.it.gam.AR <- bam(
Y ~
c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr") +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord), bs = "cr",
data = max.it,
method = "fREML",
rho = r1,
AR.start = max.it$start.event
)
summary(max.it.gam.AR)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord,
## bs = "cr") + s(X, by = vowel.ord)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.3195 0.2587 -24.428 < 2e-16 ***
## c2phonation.ordvoiceless 1.3261 0.2367 5.602 2.22e-08 ***
## c2place.ordvelar 2.3922 0.2416 9.902 < 2e-16 ***
## vowel.ordo 1.7830 0.2891 6.168 7.39e-10 ***
## vowel.ordu -0.5921 0.3083 -1.920 0.0549 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 7.950 8.534 398.352 < 2e-16 ***
## s(X):c2phonation.ordvoiceless 2.811 3.525 7.105 4.24e-05 ***
## s(X):c2place.ordvelar 8.085 8.631 96.067 < 2e-16 ***
## s(X):vowel.ordo 5.491 6.715 13.787 < 2e-16 ***
## s(X):vowel.ordu 8.673 8.958 27.925 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.749 Deviance explained = 75.1%
## fREML = 18712 Scale est. = 43.739 n = 5721
acf_plot(resid(max.it.gam.AR), split_by=list(max.it$rec.date))

acf_resid(max.it.gam.AR, split_pred = "AR.start")

plot_smooth(max.it.gam.AR, view = "X",
plot_all= "c2phonation.ord", rug = FALSE)
## Summary:
## * c2phonation.ord : factor; set to the value(s): voiced, voiceless.
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 30 values ranging from -44.040500 to 61.699900.

plot_diff(max.it.gam.AR, view = "X",
comp = list(c2phonation.ord = c("voiceless", "voiced")))
## Summary:
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 100 values ranging from -44.040500 to 61.699900.

##
## X window(s) of significant difference(s):
## -44.040500 - -7.725615
GAMMs at closure
closure.it$c2phonation.ord <- as.ordered(closure.it$c2phonation)
closure.it$c2place.ord <- as.ordered(closure.it$c2place)
closure.it$vowel.ord <- as.ordered(closure.it$vowel)
contrasts(closure.it$c2phonation.ord) <- "contr.treatment"
contrasts(closure.it$c2place.ord) <- "contr.treatment"
contrasts(closure.it$vowel.ord) <- "contr.treatment"
closure.it <- closure.it %>%
na.omit()
closure.it.gam <- bam(
Y ~
c2phonation.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr"),
data = closure.it,
method = "ML"
)
summary(closure.it.gam)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord,
## bs = "cr")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.8198 0.1835 15.371 < 2e-16 ***
## c2phonation.ordvoiceless 1.0034 0.2616 3.836 0.000126 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 8.158 8.725 584.97 < 2e-16 ***
## s(X):c2phonation.ordvoiceless 2.819 3.514 7.04 4.17e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.514 Deviance explained = 51.4%
## -ML = 33952 Scale est. = 147.58 n = 8665
closure.it.gam.null <- bam(
Y ~
s(X, bs = "cr"),
data = closure.it,
method = "ML"
)
compareML(closure.it.gam, closure.it.gam.null)
## closure.it.gam: Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord,
## bs = "cr")
##
## closure.it.gam.null: Y ~ s(X, bs = "cr")
##
## Chi-square test of ML scores
## -----
## Model Score Edf Chisq Df p.value Sig.
## 1 closure.it.gam.null 33969.15 3
## 2 closure.it.gam 33951.84 6 17.307 3.000 1.470e-07 ***
##
## AIC difference: -33.44, model closure.it.gam has lower AIC.
plot_smooth(closure.it.gam, view = "X",
plot_all= "c2phonation.ord", rug = FALSE)
## Summary:
## * c2phonation.ord : factor; set to the value(s): voiced, voiceless.
## * X : numeric predictor; with 30 values ranging from -47.595900 to 64.280100.

plot_diff(closure.it.gam, view = "X",
comp = list(c2phonation.ord = c("voiceless", "voiced")))
## Summary:
## * X : numeric predictor; with 100 values ranging from -47.595900 to 64.280100.

##
## X window(s) of significant difference(s):
## -47.595900 - -15.954203
acf_plot(resid(closure.it.gam), split_by=list(closure.it$rec.date))

closure.it.gam.2 <- bam(
Y ~
c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr") +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord), bs = "cr",
data = closure.it,
method = "fREML"
)
summary(closure.it.gam.2)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord,
## bs = "cr") + s(X, by = vowel.ord)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.3530 0.2688 8.754 < 2e-16 ***
## c2phonation.ordvoiceless 1.2010 0.2437 4.929 8.43e-07 ***
## c2place.ordvelar 0.7587 0.2475 3.066 0.00218 **
## vowel.ordo 1.9912 0.3036 6.559 5.72e-11 ***
## vowel.ordu -1.9065 0.3386 -5.630 1.86e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 6.801 7.659 329.07 < 2e-16 ***
## s(X):c2phonation.ordvoiceless 3.186 3.941 12.08 1.33e-09 ***
## s(X):c2place.ordvelar 7.512 8.274 75.43 < 2e-16 ***
## s(X):vowel.ordo 6.598 7.719 20.34 < 2e-16 ***
## s(X):vowel.ordu 8.376 8.862 25.00 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.58 Deviance explained = 58.2%
## fREML = 33340 Scale est. = 127.38 n = 8665
closure.it.gam.null <- bam(
Y ~
c2phonation.ord +
c2place.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr") +
s(X, by = c2place.ord, bs = "cr"),
data = closure.it,
method = "fREML"
)
compareML(closure.it.gam.2, closure.it.gam.null)
## closure.it.gam.2: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord,
## bs = "cr") + s(X, by = vowel.ord)
##
## closure.it.gam.null: Y ~ c2phonation.ord + c2place.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord,
## bs = "cr") + s(X, by = c2place.ord, bs = "cr")
##
## Chi-square test of fREML scores
## -----
## Model Score Edf Chisq Df p.value Sig.
## 1 closure.it.gam.null 33675.81 9
## 2 closure.it.gam.2 33339.51 15 336.302 6.000 < 2e-16 ***
##
## AIC difference: -692.69, model closure.it.gam.2 has lower AIC.
plot_smooth(closure.it.gam.2, view = "X",
plot_all= "c2phonation.ord", rug = FALSE)
## Summary:
## * c2phonation.ord : factor; set to the value(s): voiced, voiceless.
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 30 values ranging from -47.595900 to 64.280100.

plot_diff(closure.it.gam.2, view = "X",
comp = list(c2phonation.ord = c("voiceless", "voiced")))
## Summary:
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 100 values ranging from -47.595900 to 64.280100.

##
## X window(s) of significant difference(s):
## -47.595900 - -13.694082
closure.it.gamm <- bam(
Y ~
c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr") +
s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
# s(X, speaker, bs="fs", xt="cr", m=1, k=5) +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord), bs = "cr",
data = closure.it,
method = "fREML"
)
## Warning in gam.side(sm, X, tol = .Machine$double.eps^0.5): model has
## repeated 1-d smooths of same variable.
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") +
## s(X, by = vowel.ord)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.8955 1.0209 1.857 0.0634 .
## c2phonation.ordvoiceless 1.4241 0.9217 1.545 0.1224
## c2place.ordvelar 1.1709 0.9258 1.265 0.2060
## vowel.ordo 2.4567 1.1319 2.170 0.0300 *
## vowel.ordu -1.8838 1.1323 -1.664 0.0962 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 7.340 8.094 216.753 < 2e-16 ***
## s(X):c2phonation.ordvoiceless 3.641 4.424 5.772 6.97e-05 ***
## s(X,rec.date) 535.105 1255.000 7.295 < 2e-16 ***
## s(X):c2place.ordvelar 8.082 8.631 83.016 < 2e-16 ***
## s(X):vowel.ordo 7.253 8.233 13.469 < 2e-16 ***
## s(X):vowel.ordu 8.327 8.826 32.646 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.796 Deviance explained = 81%
## fREML = 30882 Scale est. = 61.807 n = 8665
closure.it.gam.null <- bam(
Y ~
c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr") +
# s(X, rec.date, bs="fs", xt="cr", m=1, k=5) +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord), bs = "cr",
data = closure.it,
method = "fREML"
)
compareML(closure.it.gamm, closure.it.gam.null)
## closure.it.gamm: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, rec.date, bs = "fs",
## xt = "cr", m = 1, k = 5) + s(X, by = c2place.ord, bs = "cr") +
## s(X, by = vowel.ord)
##
## closure.it.gam.null: Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord,
## bs = "cr") + s(X, by = vowel.ord)
##
## Chi-square test of fREML scores
## -----
## Model Score Edf Chisq Df p.value Sig.
## 1 closure.it.gam.null 33339.51 15
## 2 closure.it.gamm 30882.43 18 2457.080 3.000 < 2e-16 ***
##
## AIC difference: -5746.74, model closure.it.gamm has lower AIC.
plot_smooth(closure.it.gamm, view = "X",
plot_all= "c2phonation.ord", rug = FALSE)
## Summary:
## * c2phonation.ord : factor; set to the value(s): voiced, voiceless.
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 30 values ranging from -47.595900 to 64.280100.
## * rec.date : factor; set to the value(s): 12/12/2016 14:51:50.

plot_diff(closure.it.gamm, view = "X",
comp = list(c2phonation.ord = c("voiceless", "voiced")))
## Summary:
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 100 values ranging from -47.595900 to 64.280100.
## * rec.date : factor; set to the value(s): 12/12/2016 14:51:50.

##
## X window(s) of significant difference(s):
## -47.595900 - -19.344385
acf_plot(resid(max.it.gamm), split_by=list(max.it$rec.date))

Autoregressive correction
max.it <- arrange(max.it, rec.date, fan) %>%
mutate(rec.date = as.character(rec.date))
previous <- ""
for (i in 1:nrow(max.it)) {
if (max.it$rec.date[i] != previous) {
max.it$start.event[i] <- TRUE
previous <- max.it$rec.date[i]
} else {
max.it$start.event[i] <- FALSE
previous <- max.it$rec.date[i]
}
}
r1 <- start_value_rho(max.it.gam.2)
max.it.gam.AR <- bam(
Y ~
c2phonation.ord +
c2place.ord +
vowel.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr") +
s(X, by = c2place.ord, bs = "cr") +
s(X, by = vowel.ord), bs = "cr",
data = max.it,
method = "fREML",
rho = r1,
AR.start = max.it$start.event
)
summary(max.it.gam.AR)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + c2place.ord + vowel.ord + s(X, bs = "cr") +
## s(X, by = c2phonation.ord, bs = "cr") + s(X, by = c2place.ord,
## bs = "cr") + s(X, by = vowel.ord)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.3195 0.2587 -24.428 < 2e-16 ***
## c2phonation.ordvoiceless 1.3261 0.2367 5.602 2.22e-08 ***
## c2place.ordvelar 2.3922 0.2416 9.902 < 2e-16 ***
## vowel.ordo 1.7830 0.2891 6.168 7.39e-10 ***
## vowel.ordu -0.5921 0.3083 -1.920 0.0549 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 7.950 8.534 398.352 < 2e-16 ***
## s(X):c2phonation.ordvoiceless 2.811 3.525 7.105 4.24e-05 ***
## s(X):c2place.ordvelar 8.085 8.631 96.067 < 2e-16 ***
## s(X):vowel.ordo 5.491 6.715 13.787 < 2e-16 ***
## s(X):vowel.ordu 8.673 8.958 27.925 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.749 Deviance explained = 75.1%
## fREML = 18712 Scale est. = 43.739 n = 5721
acf_plot(resid(max.it.gam.AR), split_by=list(max.it$rec.date))

acf_resid(max.it.gam.AR, split_pred = "AR.start")

plot_smooth(max.it.gam.AR, view = "X",
plot_all= "c2phonation.ord", rug = FALSE)
## Summary:
## * c2phonation.ord : factor; set to the value(s): voiced, voiceless.
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 30 values ranging from -44.040500 to 61.699900.

plot_diff(max.it.gam.AR, view = "X",
comp = list(c2phonation.ord = c("voiceless", "voiced")))
## Summary:
## * c2place.ord : factor; set to the value(s): coronal.
## * vowel.ord : factor; set to the value(s): a.
## * X : numeric predictor; with 100 values ranging from -44.040500 to 61.699900.

##
## X window(s) of significant difference(s):
## -44.040500 - -7.725615
GAMMs at peak 1
splines.it$c2phonation.ord <- as.ordered(splines.it$c2phonation)
contrasts(splines.it$c2phonation.ord) <- "contr.treatment"
peak.it <- filter(splines.it, label %in% c("peak1_TT", "peak1_TD"))
peak.it.gam <- bam(
Y ~
c2phonation.ord +
s(X, bs = "cr") +
s(X, by = c2phonation.ord, bs = "cr"),
data = peak.it,
method = "ML"
)
summary(peak.it.gam)
##
## Family: gaussian
## Link function: identity
##
## Formula:
## Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord,
## bs = "cr")
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.8024 0.1809 -32.083 <2e-16 ***
## c2phonation.ordvoiceless 0.5019 0.2541 1.976 0.0482 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(X) 8.691 8.958 562.990 <2e-16 ***
## s(X):c2phonation.ordvoiceless 1.695 2.137 0.932 0.425
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.59 Deviance explained = 59%
## -ML = 22127 Scale est. = 96.198 n = 5971
peak.it.gam.null <- bam(
Y ~
s(X, bs = "cr"),
data = peak.it,
method = "ML"
)
compareML(peak.it.gam, peak.it.gam.null)
## peak.it.gam: Y ~ c2phonation.ord + s(X, bs = "cr") + s(X, by = c2phonation.ord,
## bs = "cr")
##
## peak.it.gam.null: Y ~ s(X, bs = "cr")
##
## Chi-square test of ML scores
## -----
## Model Score Edf Chisq Df p.value Sig.
## 1 peak.it.gam.null 22129.00 3
## 2 peak.it.gam 22126.83 6 2.171 3.000 0.227
##
## AIC difference: -0.63, model peak.it.gam has lower AIC.
## Warning in compareML(peak.it.gam, peak.it.gam.null): Only small difference in ML...
plot_smooth(peak.it.gam, view = "X",
plot_all= "c2phonation.ord", rug = FALSE)
## Summary:
## * c2phonation.ord : factor; set to the value(s): voiced, voiceless.
## * X : numeric predictor; with 30 values ranging from -48.091200 to 64.315100.

plot_diff(peak.it.gam, view = "X",
comp = list(c2phonation.ord = c("voiceless", "voiced")))
## Summary:
## * X : numeric predictor; with 100 values ranging from -48.091200 to 64.315100.

##
## Difference is not significant.